clear all;
clc;

% ------------------ Import Thresholds ------------------------------------

cd('p1_thresholds');

epith_pre = 'neg_epi_bp_dr_act_9um_';
subth_pre = 'neg_sub_bp_dr_act_9um_';

tissname = 'vert';

epith_con = load([epith_pre,'control','.txt']);
epith_low = load([epith_pre,tissname,'low','.txt']);
epith_upp = load([epith_pre,tissname,'high','.txt']);

subth_con = load([subth_pre,'control','.txt']);
subth_low = load([subth_pre,tissname,'low','.txt']);
subth_upp = load([subth_pre,tissname,'high','.txt']);

cd ..;

% ------------------ Import Unit Currents ---------------------------------

cd('p1_potentials');

epicur_pre = 'epi_bp_current_9um_';
subcur_pre = 'sub_bp_current_9um_';

epi_icon = load([epicur_pre,'control','.txt']);
epi_ilow = load([epicur_pre,tissname,'low','.txt']);
epi_iupp = load([epicur_pre,tissname,'high','.txt']);

sub_icon = load([subcur_pre,'control','.txt']);
sub_ilow = load([subcur_pre,tissname,'low','.txt']);
sub_iupp = load([subcur_pre,tissname,'high','.txt']);

cd ..;

% ------------------ Processing Data --------------------------------------

box_groups = {'0.015','0.03','0.045'};

e1 = abs( epith_low(:,1) )/(2/epi_ilow/1e3);
e2 = abs( epith_con(:,1) )/(2/epi_icon/1e3); 
e3 = abs( epith_upp(:,1) )/(2/epi_iupp/1e3);
epidata = [e1,e2,e3];

s1 = abs( subth_low(:,1) )/(2/sub_ilow/1e3);
s2 = abs( subth_con(:,1) )/(2/sub_icon/1e3); 
s3 = abs( subth_upp(:,1) )/(2/sub_iupp/1e3);
subdata = [s1,s2,s3];

isens_expsub = [0.20,0.35,0.70,0.55,0.80];
ipain_expsub = [0.50,0.50,1.50,1.40,3.10];

isens_expepi = [1.30,19.0,7.00,5.00,9.00];
ipain_expepi = [2.20,24.5,19.0,5.00,19.0];    

nogroup = length(box_groups);

x = 0:1:(nogroup+1);
isens_epiavg = median(isens_expepi)*ones(size(x));
isens_subavg = median(isens_expsub)*ones(size(x));
ipain_epiavg = median(ipain_expepi)*ones(size(x));
ipain_subavg = median(ipain_expsub)*ones(size(x));

display('50 % lower');
mean(abs(epidata(:,2)-epidata(:,1))./epidata(:,2))
(epi_icon-epi_ilow)/epi_icon

display('50 % higher');
mean(abs(epidata(:,2)-epidata(:,3))./epidata(:,2))
(epi_icon-epi_iupp)/epi_icon

% ------------------ Plotting Data ----------------------------------------

figure;
subplot(1,2,1);
hold on;
plot(x,isens_epiavg,'k--','LineWidth',4);
plot(x,ipain_epiavg,'k-','LineWidth',4);
h1 = boxplot(epidata,'labels',box_groups);
hold off;
title('Epidural','FontSize',30);
xlabel({'';'\sigma_{dura} (S/m)'},'FontSize',30);
ylabel('Stimulation Threshold Current (mA)','FontSize',30);
set(gca,'FontSize',26,'YTick',[0.1,1,10],'YScale','Log');
set(gca,'xtick',1:size(box_groups,2),'xticklabel',box_groups,'FontSize',26);
set(gca,'LineWidth',2);
set(h1,'LineWidth',3);
ylim([0.09,20]);

subplot(1,2,2);
hold on;
plot(x,isens_subavg,'k--','LineWidth',4);
plot(x,ipain_subavg,'k-','LineWidth',4);
h2 = boxplot(subdata,'labels',box_groups);
hold off;
title('Intradural','FontSize',30);
xlabel({'';'\sigma_{dura} (S/m)'},'FontSize',30);
set(gca,'FontSize',26,'YTick',[0.1,1,10],'YScale','Log');
set(gca,'xtick',1:size(box_groups,2),'xticklabel',box_groups,'FontSize',26);
set(gca,'LineWidth',2);
set(h2,'LineWidth',3);
legend('Sensory Threshold','Pain Threshold','Location','NE');
ylim([0.09,20]);

